home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / scale_geom / scale_g.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-09-17  |  28.8 KB  |  863 lines
  1. /*
  2.  * zoom: subroutines for in-place one-pass filtered image zooming
  3.  *
  4.  * Zooms (resizes) one image into another with independent control of
  5.  * scale, translation, and filtering in x and y.  Magnifies or minifies.
  6.  * Filters with an arbitrary separable finite impulse response filter.
  7.  * (A filter h(x,y) is called separable if h(x,y)=f(x)*g(y)).
  8.  * See the filt package regarding filtering.
  9.  * The code is graphics device independent, using only generic scanline
  10.  * i/o operations.
  11.  *
  12.  * The program makes one pass over the image using a moving window of
  13.  * scanlines, so memory usage is proportional to imagewidth*filterheight,
  14.  * not imagewidth*imageheight, as you'd get if the entire image were
  15.  * buffered.  The separability of the filter is also exploited, to get
  16.  * a cpu time proportional to npixels*(filterwidth+filterheight),
  17.  * not npixels*filterwidth*filterheight as with direct 2-D convolution.
  18.  *
  19.  * The source and destination images can be identical, with overlapping
  20.  * windows, in which case the algorithm finds the magic split point and
  21.  * uses the appropriate scanning directions to avoid feedback (recursion)
  22.  * in the filtering.
  23.  *
  24.  * Paul Heckbert    12 Sept 1988
  25.  * UC Berkeley
  26.  * ph@miro.berkeley.edu
  27.  *
  28.  * Copyright (c) 1989  Paul S. Heckbert
  29.  * This source may be used for peaceful, nonprofit purposes only, unless
  30.  * under licence from the author. This notice should remain in the source.
  31.  */
  32.  
  33. /*
  34.  * routines in this source file:
  35.  *
  36.  * PUBLIC:
  37.  * zoom            most common entry point; window-to-window zoom
  38.  * zoom_continuous    fancier: explicit, continuous scale&tran control
  39.  *
  40.  * SORTA PRIVATE:
  41.  * zoom_unfiltered_mono    point-sampled zooming, called by zoom_continuous
  42.  * zoom_unfiltered_rgb    point-sampled zooming, called by zoom_continuous
  43.  * zoom_filtered    filtered zooming, called by zoom_continuous
  44.  * zoom_filtered_xy    zoom, filtering x before y, called by zoom_filtered
  45.  * zoom_filtered_yx    zoom, filtering y before x, called by zoom_filtered
  46.  *
  47.  * make_weighttab    make lookup table of filter coefficients
  48.  * make_map_table    in-place tricks; order scanlines to avoid feedback
  49.  */
  50.  
  51. static char rcsid[] = "$Header: zoom.c,v 3.0 88/10/10 13:52:07 ph Locked $";
  52.  
  53. #include <math.h>
  54.  
  55. #include "simple.h"
  56. #include "pic.h"
  57. #include "filt.h"
  58. #include "scanline.h"
  59. #include "scaleg.h"
  60.  
  61. #define EPSILON 1e-7            /* error tolerance */
  62.  
  63. #define WEIGHTBITS  14            /* # bits in filter coefficients */
  64. #define FINALSHIFT  (2*WEIGHTBITS-CHANBITS) /* shift after x&y filter passes */
  65. #define WEIGHTONE  (1<<WEIGHTBITS)    /* filter weight of one */
  66.  
  67. #define INTEGER(x) (fabs((x)-floor((x)+.5)) < EPSILON)    /* is x an integer? */
  68. #define FRAC(x) fabs((x)-floor((x)+.5))        /* diff from closest integer */
  69.  
  70. typedef struct {    /* ZOOM-SPECIFIC FILTER PARAMETERS */
  71.     double scale;    /* filter scale (spacing between centers in a space) */
  72.     double supp;    /* scaled filter support radius */
  73.     int wid;        /* filter width: max number of nonzero samples */
  74. } Filtpar;
  75.  
  76. int zoom_debug = 0;    /* debug level: 0=none, 1=scanline i/o, 2=filters */
  77. int zoom_coerce = 1;    /* simplify filters if possible? */
  78. int zoom_xy = UNDEF;    /* filter x before y (1) or vice versa (0)? */
  79. int zoom_trimzeros = 1;    /* trim zeros in x filter weight tables? */
  80.  
  81. static int nzero = 0, ntot = 0, nzmax = 0;
  82.  
  83. Filt *simplify_filter();
  84.  
  85. /*----------------------------------------------------------------------
  86.  * The mapping from source coordinates to dest coordinates:
  87.  * (Notation: prefix 'a' denotes source coords, 'b' denotes destination coords)
  88.  *
  89.  *    bx = sx*ax+tx
  90.  *    by = sy*ay+ty
  91.  *
  92.  * where (ax,ay) and (bx,by) are DISCRETE source and dest coords, respectively.
  93.  *
  94.  * An important fine point on terminology: there are two kinds of pixel coords:
  95.  *    DISCRETE COORDINATES take on integer values at pixel centers
  96.  *    CONTINUOUS COORDINATES take on integer values halfway between pixels
  97.  * For example, an image with discrete coordinate range [0..511] has a
  98.  * continuous coordinate range of [0..512].  The pixel with
  99.  * discrete coords (x,y) has its center at continuous coords (x+.5,y+.5)
  100.  * and its continuous coord domain is [x..x+1] in X and [y..y+1] in Y.
  101.  * Note: discrete coords are not always stored in ints and continuous coords
  102.  * are not always stored in floats.
  103.  *
  104.  * conversion:
  105.  * if c = continuous coord and d = discrete coord, then
  106.  *    c = d+.5
  107.  *    d = floor(c)
  108.  *
  109.  * To map a discrete src interval [a0..a1] to a discrete dst interval [b0..b1]:
  110.  *
  111.  *    b-b0+.5 = bn/an*(a-a0+.5)
  112.  *  or    b = (bn/an)*a + (b0-.5)-(bn/an)*(a0-.5)
  113.  *      = scale*a+tran
  114.  *
  115.  * where a and b are the discrete source and dest coords (either x or y)
  116.  * and an and bn are the interval lengths: an=a1-a0+1, bn=b1-b0+1.
  117.  * a0, an, b0, bn are the mapping parameters used by the zoom routine below.
  118.  * In general, however, for sub-pixel scale and translate control, we allow
  119.  * any real numbers for scale and tran (although there should usually be some
  120.  * correspondence between the mapping function and the src and dst windows).
  121.  *
  122.  * We'll often want the inverse mapping, from discrete dst coord b to
  123.  * continuous src coord ac, relative to the interval origins a0 and b0:
  124.  *
  125.  *    b+b0 = s*(ac+a0-.5)+t
  126.  *  so    ac = (b+b0-s*(a0-.5)-t)/s
  127.  *       = (b+offset)/scale = MAP(b, scale, offset)
  128.  *
  129.  * The mapping fields ux and uy in the Mapping structure
  130.  * are these offsets in x & y respectively.
  131.  */
  132.  
  133. /* the mapping from discrete dest coord b to continuous source coord: */
  134. #define MAP(b, scale, offset)  ((b)+(offset))/(scale)
  135.  
  136. /*----------------------------------------------------------------------
  137.  * zoom: simplest entry point; window-to-window zoom.
  138.  * Zoom window a of apic to window b of bpic using the specified filters.
  139.  * Window a represents the pixels with x in [a->x0..a->x1]
  140.  * and y in [a->y0..a->y1], (ranges are inclusive).  Likewise for b.
  141.  * apic and bpic may be equal and a and b may overlap, no sweat.
  142.  */
  143.  
  144. zoom(apic, a, bpic, b, xfilt, yfilt)
  145. Pic *apic, *bpic;    /* source and dest pictures */
  146. Window_box *a, *b;    /* source and dest windows */
  147. Filt *xfilt, *yfilt;    /* filters for x and y */
  148. {
  149.     Mapping map;
  150.  
  151.     if (b->x0==UNDEF) {
  152.     fprintf(stderr,"zoom: undefined window for dest %s: ",
  153.         pic_get_name(bpic));
  154.     window_print("", &b);
  155.     fprintf(stderr,"\n");
  156.     exit(1);
  157.     }
  158.     window_box_set_size(a);
  159.     window_box_set_size(b);
  160.     if (a->nx<=0 || a->ny<=0) return;
  161.     map.sx = (double)b->nx/a->nx;
  162.     map.sy = (double)b->ny/a->ny;
  163.     map.tx = b->x0-.5 - map.sx*(a->x0-.5);
  164.     map.ty = b->y0-.5 - map.sy*(a->y0-.5);
  165.     zoom_continuous(apic, a, bpic, b, &map, xfilt, yfilt);
  166. }
  167.  
  168. /*
  169.  * zoom_opt: window-to-window zoom with options.
  170.  * Like zoom() but has options to make mapping square (preserve pixel shape)
  171.  * and round to the next lowest integer scale factor.
  172.  */
  173.  
  174. zoom_opt(apic, a, bpic, b, xfilt, yfilt, square, intscale)
  175. Pic *apic, *bpic;    /* source and dest pictures */
  176. Window_box *a, *b;    /* source and dest windows */
  177. Filt *xfilt, *yfilt;    /* filters for x and y */
  178. int square, intscale;    /* square mapping? integer scales? */
  179. {
  180.     Mapping map;
  181.  
  182.     if (b->x0==UNDEF) {
  183.     fprintf(stderr,"zoom_opt: undefined window for dest %s: ",
  184.         pic_get_name(bpic));
  185.     window_print("", &b);
  186.     fprintf(stderr,"\n");
  187.     exit(1);
  188.     }
  189.     window_box_set_size(a);
  190.     window_box_set_size(b);
  191.     if (a->nx<=0 || a->ny<=0) return;
  192.     map.sx = (double)b->nx/a->nx;
  193.     map.sy = (double)b->ny/a->ny;
  194.     if (square) {
  195.     if (map.sx>map.sy) {    /* use the smaller scale for both sx and sy */
  196.         map.sx = map.sy;
  197.         b->nx = ceil(a->nx*map.sx);
  198.     }
  199.     else {
  200.         map.sy = map.sx;
  201.         b->ny = ceil(a->ny*map.sy);
  202.     }
  203.     }
  204.     if (intscale) {
  205.     if (map.sx>1.) {
  206.         map.sx = floor(map.sx);
  207.         b->nx = ceil(a->nx*map.sx);
  208.     }
  209.     if (map.sy>1.) {
  210.         map.sy = floor(map.sy);
  211.         b->ny = ceil(a->ny*map.sy);
  212.     }
  213.     }
  214.     window_box_set_max(b);
  215.     map.tx = b->x0-.5 - map.sx*(a->x0-.5);
  216.     map.ty = b->y0-.5 - map.sy*(a->y0-.5);
  217.     zoom_continuous(apic, a, bpic, b, &map, xfilt, yfilt);
  218. }
  219.  
  220. /*
  221.  * zoom_continuous: zoom with explicit, continuous scale&tran control.
  222.  * Zoom according to the mapping defined in m,
  223.  * reading from window awin of apic and writing to window bwin of bpic.
  224.  * Like zoom but allows continuous, subpixel scales and translates.
  225.  * Meaning of m is explained above.  Pixels in the dest window bwin outside
  226.  * of the support of the mapped source window will be left unchanged.
  227.  * Scales in m should be positive.
  228.  */
  229.  
  230. zoom_continuous(apic, awin, bpic, bwin, m, xfilt, yfilt)
  231. Pic *apic, *bpic;
  232. Window_box *awin, *bwin;
  233. Mapping *m;
  234. Filt *xfilt, *yfilt;
  235. {
  236.     int xy;
  237.     Window_box a, b, bc, t;
  238.     Filtpar ax, ay;
  239.  
  240.     if (m->sx<=0. || m->sy<=0.) {
  241.     fprintf(stderr, "zoom_continuous: negative scales not allowed\n");
  242.     return;
  243.     }
  244.  
  245.     a = *awin;
  246.     b = *bwin;
  247.     window_box_set_size(&a);
  248.     window_box_set_size(&b);
  249.  
  250.     /*
  251.      * find scale of filter in a space (source space)
  252.      * when minifying, ascale=1/scale, but when magnifying, ascale=1
  253.      */
  254.     ax.scale = xfilt->blur*MAX(1., 1./m->sx);
  255.     ay.scale = yfilt->blur*MAX(1., 1./m->sy);
  256.     /*
  257.      * find support radius of scaled filter
  258.      * if ax.supp and ay.supp are both <=.5 then we've got point sampling.
  259.      * Point sampling is essentially a special filter whose width is fixed
  260.      * at one source pixel.
  261.      */
  262.     ax.supp = MAX(.5, ax.scale*xfilt->supp);
  263.     ay.supp = MAX(.5, ay.scale*yfilt->supp);
  264.     ax.wid = ceil(2.*ax.supp);
  265.     ay.wid = ceil(2.*ay.supp);
  266.  
  267.     /* clip source and dest windows against their respective pictures */
  268.     window_box_intersect(&a, (Window_box *)pic_get_window(apic, &t), &a);
  269.     if (pic_get_window(bpic, &t)->x0 != UNDEF)
  270.     window_box_intersect(&b, &t, &b);
  271.  
  272.     /* find valid dest window (transformed source + support margin) */
  273.     bc.x0 =  ceil(m->sx*(a.x0-ax.supp)+m->tx+EPSILON);
  274.     bc.y0 =  ceil(m->sy*(a.y0-ay.supp)+m->ty+EPSILON);
  275.     bc.x1 = floor(m->sx*(a.x1+ax.supp)+m->tx-EPSILON);
  276.     bc.y1 = floor(m->sy*(a.y1+ay.supp)+m->ty-EPSILON);
  277.     window_box_set_size(&bc);
  278.  
  279.     if (b.x0<bc.x0 || b.x1>bc.x1 || b.y0<bc.y0 || b.y1>bc.y1) {
  280.     /* requested dest window lies outside the valid dest, so clip dest */
  281.         fprintf(stderr,"\n Error: destination window not valid. \n\n");
  282.         exit(1);
  283.  
  284.       /* Note: can't clip the dest window because then the hips
  285.           header will be wrong ( Probably could re-write the header
  286.            here, but I'm not going to worry about it until a need arises -BT
  287.  
  288.     window_box_print("    clipping odst=", &b);
  289.     window_box_print(" to ", &bc);
  290.     fprintf(stderr,"\n");
  291.     window_box_intersect(&b, &bc, &b);
  292.       */
  293.     }
  294.     pic_set_window(bpic, &b);
  295.  
  296.     /* compute offsets for MAP (these will be .5 if zoom() routine was called)*/
  297.     m->ux = b.x0-m->sx*(a.x0-.5)-m->tx;
  298.     m->uy = b.y0-m->sy*(a.y0-.5)-m->ty;
  299.  
  300. #ifdef DEBUG
  301.     fprintf(stderr,"src=%s:%s", pic_get_dev(apic), pic_get_name(apic));
  302.     window_box_print("", &a);
  303.     fprintf(stderr,"\ndst=%s:%s", pic_get_dev(bpic), pic_get_name(bpic));
  304.     window_box_print("", &b);
  305.     fprintf(stderr,"\n");
  306. #endif
  307.  
  308.     if (a.nx<=0 || a.ny<=0 || b.nx<=0 || b.ny<=0) {
  309.       fprintf(stderr,"\n Error: negative window size. \n" );
  310.       return;
  311.     }
  312.  
  313.     /* check for high-level simplifications of filter */
  314.     xfilt = simplify_filter("x", zoom_coerce, xfilt, &ax, m->sx, m->ux);
  315.     yfilt = simplify_filter("y", zoom_coerce, yfilt, &ay, m->sy, m->uy);
  316.  
  317.     /*
  318.      * decide which filtering order (xy or yx) is faster for this mapping by
  319.      * counting convolution multiplies
  320.      */
  321.     xy = zoom_xy!=UNDEF ? zoom_xy :
  322.     b.nx*(a.ny*ax.wid+b.ny*ay.wid) < b.ny*(b.nx*ax.wid+a.nx*ay.wid);
  323.  
  324.     /* choose the appropriate filtering routine */
  325.     if (str_eq(xfilt->name, "point") && str_eq(yfilt->name, "point"))
  326.     zoom_unfiltered(apic, &a, bpic, &b, m);
  327.     else if (xy)
  328.     zoom_filtered_xy(apic, &a, bpic, &b, m, xfilt, &ax, yfilt, &ay);
  329.     else
  330.     zoom_filtered_yx(apic, &a, bpic, &b, m, xfilt, &ax, yfilt, &ay);
  331. }
  332.  
  333. /*
  334.  * filter_simplify: check if our discrete sampling of an arbitrary continuous
  335.  * filter, parameterized by the filter spacing (a->scale), its radius (a->supp),
  336.  * and the scale and offset of the coordinate mapping (s and u), causes the
  337.  * filter to reduce to point sampling.
  338.  *
  339.  * It reduces if support is less than 1 pixel or
  340.  * if integer scale and translation, and filter is cardinal
  341.  */
  342.  
  343. static Filt *simplify_filter(dim, coerce, filter, a, s, u)
  344. char *dim;
  345. int coerce;
  346. Filt *filter;
  347. Filtpar *a;
  348. double s, u;    /* scale and offset */
  349. {
  350.     if (coerce && (a->supp<=.5 ||
  351.     filter->cardinal && INTEGER(1./a->scale) && INTEGER(1./(s*a->scale))
  352.     && INTEGER((u/s-.5)/a->scale))) {
  353.         if (!str_eq(filter->name, "point"))
  354.         fprintf(stderr,"coercing %sfilter=%s to point\n", dim, filter->name);
  355.         filter = filt_find("point");
  356.         a->scale = 1.;
  357.         a->supp = .5;
  358.         a->wid = 1;
  359.     }
  360.     return filter;
  361. }
  362.  
  363. /*----------------------------------------------------------------------*/
  364.  
  365. /* zoom_unfiltered: do unfiltered zoom (point sampling) */
  366.  
  367. zoom_unfiltered(apic, a, bpic, b, m)
  368. Pic *apic, *bpic;
  369. Window_box *a, *b;
  370. Mapping *m;
  371. {
  372.     int overlap;
  373.  
  374.     /* do source and dest windows overlap? */
  375.     overlap = apic==bpic && window_box_overlap(a, b);
  376.  
  377. #ifdef DEBUG
  378.     fprintf(stderr,"-map (sx,sy,tx,ty)= %.2g %.2g %.2g %.2g \n",
  379.      m->sx, m->sy, m->tx, m->ty);
  380. #endif
  381.     /*
  382.      * note: We have only x-resample before y-resample versions of the
  383.      * unfiltered zoom case; we could optimize further by creating
  384.      * y-resample before x-resample versions.
  385.      */
  386.  
  387.     zoom_unfiltered_mono(apic, a, bpic, b, m, overlap);
  388. }
  389.  
  390. /* zoom_unfiltered_mono: monochrome unfiltered zoom */
  391.  
  392. zoom_unfiltered_mono(apic, a, bpic, b, m, overlap)
  393. Pic *apic, *bpic;
  394. Window_box *a, *b;
  395. Mapping *m;
  396. int overlap;
  397. {
  398.     int byi, by, ay, ayold, bx, ax;
  399.     Pixel1 *abuf, *bbuf, *bp;    /* source and dest scanline buffers */
  400.     Pixel1 **xmap, **xp;    /* mapping from abuf to bbuf (optimization) */
  401.     short *ymap;        /* order of dst y coords that avoids feedback */
  402.     char *linewritten;        /* has scanline y been written? (debugging) */
  403.  
  404.     ALLOC(abuf, Pixel1, a->nx);
  405.     ALLOC(bbuf, Pixel1, b->nx);
  406.     ALLOC(xmap, Pixel1 *, b->nx);
  407.     ALLOC(ymap, short, b->ny);
  408.     ALLOC_ZERO(linewritten, char, MAX(a->y1, b->y1)+1);
  409.  
  410.     /* if overlapping src & dst, find magic y ordering that avoids feedback */
  411.     make_map_table(m->sy, m->ty, .5, a->y0, b->y0, b->ny, overlap, ymap);
  412.  
  413.     for (bx=0; bx<b->nx; bx++) {
  414.     ax = MAP(bx, m->sx, m->ux);
  415.     xmap[bx] = &abuf[ax];
  416.     }
  417.  
  418.     ayold = -1;        /* impossible value for ay */
  419.     for (byi=0; byi<b->ny; byi++) {
  420.     by = ymap[byi];
  421.     ay = MAP(by, m->sy, m->uy);
  422.     /* scan line a.y0+ay goes to b.y0+by */
  423.     if (ay!=ayold) {        /* new scan line; read it in */
  424.         ayold = ay;
  425.         if (overlap && linewritten[a->y0+ay])
  426.         fprintf(stderr, "FEEDBACK ON LINE %d\n", a->y0+ay);
  427.         /* read scan line into abuf */
  428.         pic_read_row(apic, a->y0+ay, a->x0, a->nx, abuf);
  429.         /* resample in x */
  430.         for (bp=bbuf, xp=xmap, bx=0; bx<b->nx; bx++)
  431.         *bp++ = **xp++;
  432.     }
  433.     pic_write_row(bpic, b->y0+by, b->x0, b->nx, bbuf);
  434.     linewritten[b->y0+by] = 1;
  435.     }
  436.     free(abuf);
  437.     free(bbuf);
  438.     free(xmap);
  439.     free(ymap);
  440. }
  441. /*----------------------------------------------------------------------*/
  442.  
  443. /*
  444.  * zoom_filtered_xy: filtered zoom, xfilt before yfilt
  445.  *
  446.  * note: when calling make_weighttab, we can trim leading and trailing
  447.  * zeros from the x weight buffers as an optimization,
  448.  * but not for y weight buffers since the split formula is anticipating
  449.  * a constant amount of buffering of source scanlines;
  450.  * trimming zeros in yweight could cause feedback.
  451.  */
  452.  
  453.  
  454. zoom_filtered_xy(apic, a, bpic, b, m, xfilt, ax, yfilt, ay)
  455. Pic *apic, *bpic;
  456. Window_box *a, *b;
  457. Mapping *m;
  458. Filt *xfilt, *yfilt;    /* x and y filters */
  459. Filtpar *ax, *ay;    /* extra x and y filter parameters */
  460. {
  461.     int ayf, bx, byi, by, overlap, nchan;
  462.  
  463.                 /*PIXELTYPE NBITS  RAISON D'ETRE */
  464.     Scanline abbuf;        /*   PIXEL1   8  src&dst scanline bufs */
  465.     Scanline accum;        /*   PIXEL4  28  accumulator buf for yfilt */
  466.     Scanline *linebuf, *tbuf;    /*   PIXEL2  14  circular buf of active lines */
  467.  
  468.     short *ymap;        /* order of dst y coords that avoids feedback */
  469.     Weighttab *xweights;    /* sampled filter at each dst x pos; for xfilt*/
  470.     short *xweightbuf, *xwp;    /* big block of memory addressed by xweights */
  471.     Weighttab yweight;        /* a single sampled filter for current y pos */
  472.     char *linewritten;        /* has scanline y been written? (debugging) */
  473.  
  474.     nchan = pic_get_nchan(apic)==1 ? PIXEL_MONO : PIXEL_RGB;
  475.     scanline_alloc(&abbuf, PIXEL1|nchan, MAX(a->nx, b->nx));
  476.     scanline_alloc(&accum, PIXEL4|nchan, b->nx);
  477.     ALLOC(linebuf, Scanline, ay->wid);
  478.     /* construct circular buffer of ay->wid intermediate scanlines */
  479.     for (ayf=0; ayf<ay->wid; ayf++) {
  480.     scanline_alloc(&linebuf[ayf], PIXEL2|nchan, b->nx);
  481.     linebuf[ayf].y = -1;        /* mark scanline as unread */
  482.     }
  483.  
  484.     ALLOC(ymap, short, b->ny);
  485.     ALLOC(xweights, Weighttab, b->nx);
  486.     ALLOC(xweightbuf, short, b->nx*ax->wid);
  487.     ALLOC(yweight.weight, short, ay->wid);
  488.     ALLOC_ZERO(linewritten, char, MAX(a->y1, b->y1)+1);
  489.  
  490.     /* do source and dest windows overlap? */
  491.     overlap = apic==bpic && window_box_overlap(a, b);
  492.  
  493. #ifdef DEBUG
  494.     fprintf(stderr,"-xy -map %g %g %g %g\n", m->sx, m->sy, m->tx, m->ty);
  495.     fprintf(stderr,"X: filt=%s support=%g,  scale=%g scaledsupport=%g width=%d\n",
  496.     xfilt->name, xfilt->supp, ax->scale, ax->supp, ax->wid);
  497.     fprintf(stderr,"Y: filt=%s support=%g,  scale=%g scaledsupport=%g width=%d\n",
  498.     yfilt->name, yfilt->supp, ay->scale, ay->supp, ay->wid);
  499. #endif
  500.  
  501.     /*
  502.      * prepare a weighttab (a sampled filter for source pixels) for
  503.      * each dest x position
  504.      */
  505.     for (xwp=xweightbuf, bx=0; bx<b->nx; bx++, xwp+=ax->wid) {
  506.     xweights[bx].weight = xwp;
  507.     make_weighttab(b->x0+bx, MAP(bx, m->sx, m->ux),
  508.         xfilt, ax, a->nx, zoom_trimzeros, &xweights[bx]);
  509.     }
  510.  
  511.     /* if overlapping src & dst, find magic y ordering that avoids feedback */
  512.     make_map_table(m->sy, m->ty, ay->supp, a->y0, b->y0, b->ny, overlap,
  513.     ymap);
  514.  
  515.     for (byi=0; byi<b->ny; byi++) {    /* loop over dest scanlines */
  516.     by = ymap[byi];
  517.     if (zoom_debug) fprintf(stderr,"by=%d: ", b->y0+by);
  518.     /* prepare a weighttab for dest y position by */
  519.     make_weighttab(b->y0+by, MAP(by, m->sy, m->uy),
  520.         yfilt, ay, a->ny, 0, &yweight);
  521.     if (zoom_debug) fprintf(stderr,"ay=%d-%d, reading: ",
  522.         a->y0+yweight.i0, a->y0+yweight.i1-1);
  523.     scanline_zero(&accum);
  524.  
  525.     /* loop over source scanlines that influence this dest scanline */
  526.     for (ayf=yweight.i0; ayf<yweight.i1; ayf++) {
  527.         tbuf = &linebuf[ayf % ay->wid];
  528.         if (tbuf->y != ayf) {    /* scanline needs to be read in */
  529.         if (zoom_debug) fprintf(stderr,"%d ", a->y0+ayf);
  530.         if (overlap && linewritten[a->y0+ayf])
  531.             fprintf(stderr, "FEEDBACK ON LINE %d\n", a->y0+ayf);
  532.         tbuf->y = ayf;
  533.         /* read source scanline into abbuf */
  534.         abbuf.len = a->nx;
  535.         scanline_read(apic, a->x0, a->y0+ayf, &abbuf);
  536.         /* and filter it into the appropriate line of linebuf (xfilt) */
  537.         scanline_filter(CHANBITS, xweights, &abbuf, tbuf);
  538.         }
  539.         /* add weighted tbuf into accum (these do yfilt) */
  540.         scanline_accum(yweight.weight[ayf-yweight.i0], tbuf, &accum);
  541.     }
  542.  
  543.     /* shift accum down into abbuf and write out dest scanline */
  544.     abbuf.len = b->nx;
  545.     scanline_shift(FINALSHIFT, &accum, &abbuf);
  546.     scanline_write(bpic, b->x0, b->y0+by, &abbuf);
  547.     linewritten[b->y0+by] = 1;
  548.     if (zoom_debug) fprintf(stderr,"\n");
  549.     }
  550.  
  551.     scanline_free(&abbuf);
  552.     scanline_free(&accum);
  553.     for (ayf=0; ayf<ay->wid; ayf++)
  554.     scanline_free(&linebuf[ayf]);
  555.     free(ymap);
  556.     free(linebuf);
  557.     free(xweights);
  558.     free(xweightbuf);
  559.     free(yweight.weight);
  560.     free(linewritten);
  561.     statistics();
  562. }
  563.  
  564. /* zoom_filtered_yx: filtered zoom, yfilt before xfilt */
  565.  
  566. zoom_filtered_yx(apic, a, bpic, b, m, xfilt, ax, yfilt, ay)
  567. Pic *apic, *bpic;
  568. Window_box *a, *b;
  569. Mapping *m;
  570. Filt *xfilt, *yfilt;    /* x and y filters */
  571. Filtpar *ax, *ay;    /* extra x and y filter parameters */
  572. {
  573.     int ayf, bx, byi, by, overlap, nchan;
  574.  
  575.                 /*PIXELTYPE NBITS  RAISON D'ETRE */
  576.     Scanline bbuf;        /*   PIXEL1   8  dst scanline buf */
  577.     Scanline accum;        /*   PIXEL4  22  accumulator buf for yfilt */
  578.     Scanline *linebuf, *tbuf;    /*   PIXEL1   8  circular buf of active lines */
  579.  
  580.     short *ymap;        /* order of dst y coords that avoids feedback */
  581.     Weighttab *xweights;    /* sampled filter at each dst x pos; for xfilt*/
  582.     short *xweightbuf, *xwp;    /* big block of memory addressed by xweights */
  583.     Weighttab yweight;        /* a single sampled filter for current y pos */
  584.     char *linewritten;        /* has scanline y been written? (debugging) */
  585.  
  586.     nchan = pic_get_nchan(apic)==1 ? PIXEL_MONO : PIXEL_RGB;
  587.     scanline_alloc(&bbuf, PIXEL1|nchan, b->nx);
  588.     scanline_alloc(&accum, PIXEL4|nchan, a->nx);
  589.     ALLOC(linebuf, Scanline, ay->wid);
  590.     /* construct circular buffer of ay->wid intermediate scanlines */
  591.     for (ayf=0; ayf<ay->wid; ayf++) {
  592.     scanline_alloc(&linebuf[ayf], PIXEL1|nchan, a->nx);
  593.     linebuf[ayf].y = -1;        /* mark scanline as unread */
  594.     }
  595.  
  596.     ALLOC(ymap, short, b->ny);
  597.     ALLOC(xweights, Weighttab, b->nx);
  598.     ALLOC(xweightbuf, short, b->nx*ax->wid);
  599.     ALLOC(yweight.weight, short, ay->wid);
  600.     ALLOC_ZERO(linewritten, char, MAX(a->y1, b->y1)+1);
  601.  
  602.     /* do source and dest windows overlap? */
  603.     overlap = apic==bpic && window_box_overlap(a, b);
  604.  
  605. #ifdef DEBUG
  606.     fprintf(stderr,"-yx -map %g %g %g %g\n", m->sx, m->sy, m->tx, m->ty);
  607.     fprintf(stderr,"X: filt=%s supp=%g  scale=%g scaledsupp=%g wid=%d\n",
  608.     xfilt->name, xfilt->supp, ax->scale, ax->supp, ax->wid);
  609.     fprintf(stderr,"Y: filt=%s supp=%g  scale=%g scaledsupp=%g wid=%d\n",
  610.     yfilt->name, yfilt->supp, ay->scale, ay->supp, ay->wid);
  611. #endif
  612.  
  613.     /*
  614.      * prepare a weighttab (a sampled filter for source pixels) for
  615.      * each dest x position
  616.      */
  617.     for (xwp=xweightbuf, bx=0; bx<b->nx; bx++, xwp+=ax->wid) {
  618.     xweights[bx].weight = xwp;
  619.     make_weighttab(b->x0+bx, MAP(bx, m->sx, m->ux),
  620.         xfilt, ax, a->nx, zoom_trimzeros, &xweights[bx]);
  621.     }
  622.  
  623.     /* if overlapping src & dst, find magic y ordering that avoids feedback */
  624.     make_map_table(m->sy, m->ty, ay->supp, a->y0, b->y0, b->ny, overlap,
  625.     ymap);
  626.  
  627.     for (byi=0; byi<b->ny; byi++) {     /* loop over dest scanlines */
  628.     by = ymap[byi];
  629.     if (zoom_debug) fprintf(stderr,"by=%d: ", b->y0+by);
  630.     /* prepare a weighttab for dest y position by */
  631.     make_weighttab(b->y0+by, MAP(by, m->sy, m->uy),
  632.         yfilt, ay, a->ny, 0, &yweight);
  633.     if (zoom_debug) fprintf(stderr,"ay=%d-%d, reading: ",
  634.         a->y0+yweight.i0, a->y0+yweight.i1-1);
  635.     scanline_zero(&accum);
  636.  
  637.     /* loop over source scanlines that influence this dest scanline */
  638.     for (ayf=yweight.i0; ayf<yweight.i1; ayf++) {
  639.         tbuf = &linebuf[ayf % ay->wid];
  640.         if (tbuf->y != ayf) {    /* scanline needs to be read in */
  641.         if (zoom_debug) fprintf(stderr,"%d ", a->y0+ayf);
  642.         if (overlap && linewritten[a->y0+ayf])
  643.             fprintf(stderr, "FEEDBACK ON LINE %d\n", a->y0+ayf);
  644.         tbuf->y = ayf;
  645.         /* read source scanline into linebuf */
  646.         scanline_read(apic, a->x0, a->y0+ayf, tbuf);
  647.         }
  648.         /* add weighted tbuf into accum (these do yfilt) */
  649.         scanline_accum(yweight.weight[ayf-yweight.i0], tbuf, &accum);
  650.     }
  651.  
  652.     /* and filter it into the appropriate line of linebuf (xfilt) */
  653.     scanline_filter(FINALSHIFT, xweights, &accum, &bbuf);
  654.     /* and write out dest scanline in bbuf */
  655.     scanline_write(bpic, b->x0, b->y0+by, &bbuf);
  656.     linewritten[b->y0+by] = 1;
  657.     if (zoom_debug) fprintf(stderr,"\n");
  658.     }
  659.  
  660.     scanline_free(&bbuf);
  661.     scanline_free(&accum);
  662.     for (ayf=0; ayf<ay->wid; ayf++)
  663.     scanline_free(&linebuf[ayf]);
  664.     free(ymap);
  665.     free(linebuf);
  666.     free(xweights);
  667.     free(xweightbuf);
  668.     free(yweight.weight);
  669.     free(linewritten);
  670.     statistics();
  671. }
  672.  
  673. /*
  674.  * make_weighttab: sample the continuous filter, scaled by ap->scale and
  675.  * positioned at continuous source coordinate cen, for source coordinates in
  676.  * the range [0..len-1], writing the weights into wtab.
  677.  * Scale the weights so they sum to WEIGHTONE, and trim leading and trailing
  678.  * zeros if trimzeros!=0.
  679.  * b is the dest coordinate (for diagnostics).
  680.  */
  681.  
  682. static make_weighttab(b, cen, filter, ap, len, trimzeros, wtab)
  683. int b, len;
  684. double cen;
  685. Filt *filter;
  686. Filtpar *ap;
  687. int trimzeros;
  688. Weighttab *wtab;
  689. {
  690.     int i0, i1, i, sum, t, stillzero, lastnonzero, nz;
  691.     short *wp;
  692.     double den, sc, tr;
  693.  
  694.     /* find the source coord range of this positioned filter: [i0..i1-1] */
  695.     i0 = cen-ap->supp+.5;
  696.     i1 = cen+ap->supp+.5;
  697.     if (i0<0) i0 = 0;
  698.     if (i1>len) i1 = len;
  699.     if (i0>=i1) {
  700.     fprintf(stderr, "make_weighttab: null filter at %d\n", b);
  701.     exit(1);
  702.     }
  703.     /* the range of source samples to buffer: */
  704.     wtab->i0 = i0;
  705.     wtab->i1 = i1;
  706.  
  707.     /* find scale factor sc to normalize the filter */
  708.     for (den=0, i=i0; i<i1; i++)
  709.     den += filt_func(filter, (i+.5-cen)/ap->scale);
  710.     /* set sc so that sum of sc*func() is approximately WEIGHTONE */
  711.     sc = den==0. ? WEIGHTONE : WEIGHTONE/den;
  712.     if (zoom_debug>1) fprintf(stderr,"    b=%d cen=%g scale=%g [%d..%d) sc=%g:  ",
  713.     b, cen, ap->scale, i0, i1, sc);
  714.  
  715.     /* compute the discrete, sampled filter coefficients */
  716.     stillzero = trimzeros;
  717.     for (sum=0, wp=wtab->weight, i=i0; i<i1; i++) {
  718.  
  719.     /* evaluate the filter function: */
  720.     tr = sc*filt_func(filter, (i+.5-cen)/ap->scale);
  721.  
  722.     if (tr<MINSHORT || tr>MAXSHORT) {
  723.         fprintf(stderr, "tr=%g at %d\n", tr, b);
  724.         exit(1);
  725.     }
  726.     t = floor(tr+.5);
  727.     if (stillzero && t==0) i0++;    /* find first nonzero */
  728.     else {
  729.         stillzero = 0;
  730.         *wp++ = t;            /* add weight to table */
  731.         sum += t;
  732.         if (t!=0) lastnonzero = i;    /* find last nonzero */
  733.     }
  734.     }
  735.     ntot += wtab->i1-wtab->i0;
  736.     if (sum==0) {
  737.     nz = wtab->i1-wtab->i0;
  738.     /* fprintf(stderr, "sum=0 at %d\n", b); */
  739.     wtab->i0 = wtab->i0+wtab->i1 >> 1;
  740.     wtab->i1 = wtab->i0+1;
  741.     wtab->weight[0] = WEIGHTONE;
  742.     }
  743.     else {
  744.     if (trimzeros) {        /* skip leading and trailing zeros */
  745.         /* set wtab->i0 and ->i1 to the nonzero support of the filter */
  746.         nz = wtab->i1-wtab->i0-(lastnonzero-i0+1);
  747.         wtab->i0 = i0;
  748.         wtab->i1 = i1 = lastnonzero+1;
  749.     }
  750.     else                /* keep leading and trailing zeros */
  751.         nz = 0;
  752.     if (sum!=WEIGHTONE) {
  753.         /*
  754.          * Fudge the center slightly to make sum=WEIGHTONE exactly.
  755.          * Is this the best way to normalize a discretely sampled
  756.          * continuous filter?
  757.          */
  758.         i = cen+.5;
  759.         if (i<i0) i = i0; else if (i>=i1) i = i1-1;
  760.         t = WEIGHTONE-sum;
  761.         if (zoom_debug>1) fprintf(stderr,"[%d]+=%d ", i, t);
  762.         wtab->weight[i-i0] += t;    /* fudge center sample */
  763.     }
  764.     }
  765.     if (nz>nzmax) nzmax = nz;
  766.     nzero += nz;
  767.     if (zoom_debug>1) {
  768.     fprintf(stderr,"\t");
  769.     for (wp=wtab->weight, i=i0; i<i1; i++, wp++)
  770.         fprintf(stderr,"%5d ", *wp);
  771.     fprintf(stderr,"\n");
  772.     }
  773. }
  774.  
  775. /*
  776.  * make_map_table:
  777.  * construct a table for the mapping a = (b-t)/s for b-b0 in [0..bn-1]
  778.  * ordered so that the buffer resampling operation
  779.  *        for (bi=0; bi<bn; bi++) {
  780.  *        b = map[bi];
  781.  *        a = MAP(b, scale, offset);
  782.  *        buf[b] = buf[a];
  783.  *        }
  784.  * can work IN PLACE without feedback
  785.  *
  786.  * a and b here are source and dest coordinates, in either x or y.
  787.  * This is needed only if there is overlap between source and dest windows.
  788.  */
  789.  
  790. static make_map_table(scale, tran, asupp, a0, b0, bn, overlap, map)
  791. double scale, tran;    /* scale and translate */
  792. double asupp;        /* filter support radius in source space */
  793. int a0, b0, bn;        /* source and dest positions, dest length */
  794. int overlap;        /* do source and dest overlap? */
  795. short *map;        /* mapping table */
  796. {
  797.     int split, b, i0;
  798.     double z, u;
  799.  
  800.     /* find fixed point of mapping; let split=b-b0 at the point where a=b */
  801.  
  802.     if (overlap) {            /* source and dest windows overlap */
  803.     if (scale==1.)            /* no scale change, translation only */
  804.         /* if moving left then scan to right, and vice versa */
  805.         split = a0<b0 ? 0 : bn;
  806.     else {                /* magnify or minify */
  807.  
  808.         /* THE MAGIC SPLIT FORMULA: */
  809.  
  810.         if (scale>1.) {        /* magnify */
  811.         /*
  812.          * choose integer nearest midpoint of valid interval:
  813.          *    x < split <= x+s/(s-1)
  814.          * where x=(tran+scale*asupp)/(1-scale)-b0
  815.          */
  816.         split = floor((tran+scale*asupp-.5)/(1.-scale)-b0+1.);
  817.         }
  818.         else {            /* minify */
  819.         /*
  820.          * only one valid split pt in this case:
  821.          *    x <= split < x+1
  822.          * so we take extra care (x as above)
  823.          */
  824.         split = ceil((tran+scale*asupp)/(1.-scale)-b0);
  825.         /*
  826.          * The above formula is perfect for exact arithmetic,
  827.          * but hardware roundoff can cause split to be one too big.
  828.          * Check for roundoff by simulating precisely the calculation
  829.          * of i0 in make_weighttab.
  830.          */
  831.         u = b0-scale*(a0-.5)-tran;    /* recalculate offset */
  832.         z = MAP(split-1, scale, u);    /* cen at b=split-1 */
  833.         z = z-asupp+.5;
  834.         i0 = z;                /* i0 at b=split-1 */
  835.         if (a0+i0>=b0+split)        /* feedback at b=split-1? */
  836.             split--;            /* correct for roundoff */
  837.         }
  838.         if (split<0) split = 0;
  839.         else if (split>bn) split = bn;
  840.         fprintf(stderr,"split at y=%d\n", b0+split);
  841.     }
  842.  
  843.     if (scale>=1.) {        /* magnify: scan in toward split */
  844.         for (b=0;    b<split;  b++) *map++ = b;
  845.         for (b=bn-1; b>=split; b--) *map++ = b;
  846.     }
  847.     else {                /* minify: scan out away from split */
  848.         for (b=split-1; b>=0;  b--) *map++ = b;
  849.         for (b=split;   b<bn;  b++) *map++ = b;
  850.     }
  851.     }
  852.     else                /* no overlap; either order will work */
  853.     for (b=0; b<bn; b++) *map++ = b;    /* we opt for forward order */
  854. }
  855.  
  856. static statistics()
  857. {
  858. #   ifdef DEBUG
  859.     fprintf(stderr,"%d/%d=%.2f%% of filter coeffs were zero, nzmax=%d\n",
  860.         nzero, ntot, 100.*nzero/ntot, nzmax);
  861. #   endif
  862. }
  863.